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Abstract. - We provide a clear evidence that a two species mesoscopic Lattice Boltzmann (LB) 
model with competing short-range attractive and mid-range repulsive interactions supports emer- 
gent Herschel-Bulkley (HB) rheology, i.e. a power-law dependence of the shear-stress as a function 
of the strain rate, beyond a given yield-stress threshold. This kinetic formulation supports a seam- 
less transition from fiowing to non-flowing behaviour, through a smooth tuning of the parameters 
governing the mesoscopic interactions between the two species. The present model may become a 
valuable computational tool for the investigation of the rheology of soft-glassy materials on scales 
of experimental interest. 
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The rheology of soft flowing systems, such as emul- 
sions, foams, pastes gels, and other types of complex flu- 
ids, plays a major role in modern materials science, both 
on account of its broad range of practical applications 
and because of the challenges it poses to modern non- 
equilibrium thermodynamics [1"4 . Soft-glassy materials 
of assorted nature, emulsions, foams, pastes and granular 
materials, are known to exhibit a fairly rich and complex 
rheology. Among other signatures of complex behaviour, 
such as anomalous relaxation, dynamical arrest and reflu- 
idization, stick and slip motion, the rheology of soft-glassy 
materials is often characterized by a non-linear relation be- 
tween the applied stress and the resulting strain. A popu- 
lar expression of such non-linear behaviour is provided by 
the Herschel-Bulkley (HB) relation [SHS], a = ay + AS^^, 
where a is the applied stress, S the resulting shear (in- 
verse time) and A a material constant. The HB relation is 
characterized by a non-zero yield-stress, cry, below which 
no flow takes place, and by a scaling exponent /3 ^ 1. 
Although non-linear rheological behaviour is well docu- 
mented in several experimental studies, its microscopic 
foundations still elude a thorough theoretical understand- 
ing, thereby holding back many important applications in 
fluid mechanics, material science and biology. As for most 
complex states of matter, the experimental and theoreti- 
cal investigation of soft-glassy materials draws substantial 
benefits from the additional insights provided by computer 



simulations. Simulation methods split into two major fam- 
ilies: macroscopic/continuum and microscopic/atomistic. 
The former are computationally efhcient, but require a 
certain fore-knowledge of the basic physics in order to 
supply, upfront, constitutive equations and boundary con- 
ditions. Microscopic methods require much less coarse- 
graining, and, as a consequence, less parametric input, 
but must face with a much higher computational demand. 
A third option is offered by mesoscopic methods, which, 
as implied by their very name, work at an intermediate 
level, hopefully achieving an optimal tradeoff between the 
aforementioned two. Mesoscopic models supporting HB 
rheology have been in existence for a while in the soft 
glassy materials literature, to begin with the well-known 
model by SoUich et ai, in which the authors postulate 
a model kinetic equation for the probability P{l,E,t) of 
finding a given mesoscopic region of the flow at time t, 
with a local strain I and a local maximal yield elasticity 
E [OUllj. More recently, kinetic models for the elasto- 
plastic dynamics of jammed materials, taking the form of 
non-local Boltzmann equations for the stress distribution 
function have also been proposed [12) . In this work, we 
provide the first evidence that a mesoscopic Lattice Boltz- 
mann (LB) model with competing short-range attractive 
and mid-range repulsive interactions supports emergent 
Herschel-Bulkley (HB) rheology, i.e. a power-law depen- 
dence of the shear-stress as a function of the strain rate, 
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shift Fst/ps in the baricentric velocity in ([1]). The force 
within each species, F^^ consists of an attractive (a) com- 
ponent , acting only on the first Brillouin region (bi, index 
1 — 8 in figured]), and a repulsive (r) one, acting on both 
belts (62, index f — 24 in figured]), whereas the force be- 
tween different species {X) is short-ranged and repulsive 
(acting again on the first Brillouin region): 



F,(r, t) = F^if, t) + Fl{r, t) + F^ (r, t) 
with the general structure of the forcings given by 

ie&i,2 



(2) 



Fig. 1: The set of 25 discrete velocities, including a rest particle 
(0). The first belt (1 — 8) hosts attractive A A and BB and 
repulsive AB interactions. Both belts (1 — 24) host repulsive 
AA and BB interactions. Full details coan be found in 1141. 



beyond a given yield-stress threshold. The kinetic equa- 
tion describing the fluid rheology is not postulated on the 
basis of informed insights on the physics under inspection, 
but results instead from a lattice transcription of a basic 
Boltzmann kinetic equation, equipped with some minimal 
ingredients required to reproduce the hydrodynamics of 
non-ideal fluid mixtures [TSllMl. 
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Lattice Boltzmann with multirange interactions. 

— Our system is described by a lattice version of the 
Boltzmann kinetic equation for a multicomponent fluid 
[TSHTO] with two species {A, B): 



--[fUr,t)-ft\ps,u + rFs/ps)l 



s = A,B (1) 



where fis{f,t) is the probability density function of find- 
ing a particle of species s — A, B a.t site r and time t, 
moving along the i-th lattice direction defined by the dis- 
crete speeds Ci with i = 0, 1, 6 = 24 (see figured]). The 
left hand-side of ([T]) stands for molecular free-streaming, 
whereas the right-hand side represents the coUisional re- 
laxation towards a local equilibrium /f^^'^'' (ps , m) on a time 
scale r. The equilibrium for the s specie is a function of 
the local species density (one for each species) and of the 
baricentric velocity u: 
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ft'\ps,u)=w,Ps ( 1 + 
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where c| = J2i '^i'^fx the square of the sound speed ve- 
locity, / is the unit tensor and Wi 's are equilibrium weights 
used to enforce isotropy of the hydrodynamic equations 
|f 5| . Intermolecular forces are incorporated within the 
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-psir,t)'YwiCiPs'{f+Ci,t), s ^ s' 
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with Wi the standard weights of the two-dimensional nine- 
speed lattice, G's the strength parameters. The pseudo- 
potential ips (^7 1) 1 has been taken for both species in the 
form originally suggested by Shan & Chen |16| . namely 
^'(ps) = po{l — 6"^=/^"). The parameter po is a refer- 
ence density beyond which self-interactions become van- 
ishingly small, thereby preventing mass density collapse 
(i.e. Ps 00) due to attractive interactions. Two-belt, 
(intra-species) self-interactions are introduced to allow a 
separate control of the equation of state and surface ten- 
sion, independently. In particular, one can show that, for 
a flat A/B interface, the surface tension scales like: 
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VpA ■ VpB dy. 



where the coordinate y runs across the interface and 
Gs = G^s + ^G^g. For repulsive interactions, {Gab > 0), 
the second integral at the rhs is positive-definite, since 
\7pA ■ VpB < 0. By choosing Gg > 0, the first integral 
is negative-definite and consequently one can decrease the 
surface tension by simply increasing pQ. Full details can 
be found in |I4| . As is well known, non-trivial rheologi- 
cal behaviour has been obtained by molecular dynamics 
simulation models [2TJ[22]. A basic lesson learned from 
these models is that by taking two fluids with suitable in- 
teraction parameters (involving frustration), one is able to 
observe a phenomenology in reasonable qualitative agree- 
ment with experimental results. This suggests the possi- 
bility of formulating an equivalent model at the level of a 
suitably extended kinetic Boltzmann equation with mini- 
mal ingredients (two species plus frustration) to support 
non-linear rheology. This is exactly what characterizes 
our model. The present LB scheme embeds the universal- 
ity of the conservation laws underlying the fluid equations, 
be they ideal or interacting (non-ideal), within a compu- 
tationally efficient theoretical framework. We note that 
the Shan-Chen formulation is basically an effective one- 
body closure of the many-body Liouville equation, encod- 
ing the basic symmetries of potential energy interactions 
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within a minimal lattice formulation, i.e. a one-parameter, 
nearest-neighbor, pseudo-potential. The reason why our 
model can incorporate substantial new non-ideal physics 
without taxing computational efficiency, is again univer- 
sality: once the proper competing mechanisms are put in 
place, the specific form of the interactions is largely imma- 
terial to the large-scale behaviour of the non-ideal fluid. 
Consequently, a minimal lattice pseudo-potential is suffi- 
cient. We remark that a unique feature of the present LB 
scheme, is the capability of incorporating non-linear hy- 
drodynamics nearly " for- free" , through a simple quadratic 
dependence of the local equilibria on the local flow field. 
Thanks to this property, our model can seamlessy straddle 
across various non-trivial flow regimes (flowing/arrested) 
through a smooth change of the interaction parameters. 

Numerical Results. — The computational domain is 
a square box of size Lx L covered by x Ny = 512 x 512 
lattice sites with a uniform lattice spacing dx — I. The 
simulations, performed on latest generation Graphics Pro- 
cessing Units (GPU) ^20j, require few hours for one million 
time-steps, the typical time-span of a run. With a fixed 
set of following baseline coupling parameters 0, the ref- 
erence density po is varied between po = 0.70 and 0.90 
that corresponds to a decrease of surface tension from or- 
dinary values to an almost vanishing value for po = 0.90 
(based on the use of equation (63) in [M]). The fluid 
is initialized with pA{x,y) = 0.61(1 + 0.1 sin(/ci„j/)) and 
PBix,y) = 0.61(1 + 0.01 sinihny)) with fc„ = and 
is subject to an external periodic forcing in the x direc- 
tion of the form Fx{y) = Fq SY!i{kfy), with wavenumber 
kf = The forcing amplitude Fq is tuned in such a 
way as to produce, in standard stationary flow conditions, 
a sinusoidal Kolmogorov flow of maximum speed Uo, i.e. 
Ux{x,y) = UoSui{kfy). In a previous work [M], the sys- 
tem response was monitored using the following response 

function: R(t) ~ = where Uit) is the Fourier 

transform of the line-averaged speed along the x direction 

Ex Ux{x,y,t) 
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In the above, v is the nominal kinematic viscosity of 
both fluids and D defines the effective viscosity of the two- 
fluids system. By construction, under undisturbed flow 
conditions, R — 1, so that R <^ 1 provides a direct mea- 
sure of slowing-down through enhanced effective viscosity. 
The parameter R is thus a direct measure of the effective 
fluidity of the system. However, since we are focusing on 



^We have chosen 



7.1, G^^ = 0.045 



Negative/positive signs standing for repul- 
sion/attraction, respectively, secure that both A and B fluids are 
in the liquid phase. For all simulations we have chosen a constant 
relaxation time r = 1.0. The use of a coupling-dependent relaxation 
time has never been explored in the literature and surely deserves a 
separate study on its own. 



a non-Newtonian behaviour, it proves more informative to 
inspect first the actual space-time averaged velocity pro- 
files 

f/(2/)-^^ uAy,t)dt, r»i 

as a function of the reference density po at a given forcing 
intensity with Uq = 0.1 in computational units. From 
figure [21 a flattening of the velocity profile in the central 
region of the flow is clearly observed, for all values of po > 
0.70. This is a well-known signature of non-Newtonian 
behaviour 1231. 
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Fig. 2: The average (time and x direction) velocity profile U (y) 
for different values of po: by increasing po (i-e. decreasing the 
surface tension) the velocity profile becomes flatter and with 
a lower amplitude. In the inset, we show U{y) for po = 0.79, 
Po = 0.81 and po = 0.83 rescaled in such a way that the velocity 
gradient at ?/ = is kept fixed. The inset shows that the 
different velocity profiles cannot be superimposed by a mere 
rescaling factor. 

A typical density contour of fluid A is shown in figure [3l 
To inspect the non-Newtonian behaviour on more quan- 
titative grounds, we have measured the effective viscosity 
through the ratio of the nominal shear for a standard flow, 
to the value of the shear S{y) = dU/dy provided by the 
simulation (see inset of figure [5]). At statistical steady 
state, the momentum balance equation yields dyP^y = 
Fx{y) (derivatives along x are zero by homogeneity). Inte- 
grating along y, we obtain Pxyijj) = Fx{y')dy' , which is 
known exactly at each location y, since the right-hand-side 
is nothing but the expression of the forcing. The resulting 
shear is simply collected as the spatial derivative of the 
time averaged velocity field, i.e. Sxy{y) = dyU{y). Figure 
m shows the scatter-plot of the stress a = P^y versus the 
shear S = Sxy for each value of y. This figure carries the 
central result of this work. First, it is seen that the fluid 
starts to flow only above a critical threshold (yield-stress) 
of the order of ay ^ 1-5 10""*, which is comparable with 
the maximum applied stress ctq = puUq2it/ L ~ 2.5 x 10~^. 
Remarkably, the various data, corresponding to different 
values of the forcing, all fall within basically the same mas- 
ter curve. In the lower inset, we report the flt exposing the 
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Fig. 4: The average stress Pxy{y) as a function of the observed 
shear S(y) obtained from the numerical simulation at po = 0.83 
and different values of the forcing. In the inset, we show the 
fit of P^y by using the HB form A + B for po = 0.83 and 
po = 0.81. 



Fig. 3: A typical snapshot of the fluid A density contour. 
Blue/yellow colors code for low/high density regions. Low/high 
density regions of fluid A are filled with correspondingly 
high/low density regions of fluid B. 



exponent of the HB-like relation a = ay + B , which 
yields /3 ~ 0.25, in a reasonable good agreement with pre- 
vious models [24l[25]. The upper inset shows the same fit 
for Po = 0.81, which again yields HB behaviour, although 
with a larger exponent /? ^ 0.5. 

Since our data support HB behaviour with a surface- 
tension dependent exponent, it is worth inspecting the ef- 
fect of lowering the surface tension, and eventually taking 
it nominally below zero. To this purpose, we measure the 
time-averaged response function R for different values of 
Po- Figure [5] shows a neat divergence of the reciprocal re- 
sponse function as the condition of zero-flow (total arrest) 
is approached. Incidentally, the functional dependence of 
the time averaged response function, R(po), can be fit- 
ted reasonably well by a Vogel-Fulcher-Tammann (VFT) 
law [26H28], -R^Hpo) = A exp{-^), with A = 3.51 
and B = 0.045, although other functional forms compat- 
ible with finite-density divergence cannot be ruled out. 
For instance, the value of the maximum mean velocity 
shown in the inset of figure [5l would support a simpler 
^ {p — Pc)~^ divergence. Leaving this question to a 
future and separate investigation, here we simply observe 
that the system appears to come to a complete arrest as 
the surface tension is sent to smaller and smaller values 
(the nominal zero-point is at po 0.87). Finally, we point 
out that the system can also be taken to virtually negative 
surface tensions, in which case lamellar-like configurations 
are observed. However, the physical viability/reliability of 
the present model in this parameter regime still needs to 
be assessed. 



Conclusions and Outlook. — Summarizing, we have 
provided the first evidence of emergent Herschel-Bulkley 
(HB) rheology from a "first principle" lattice kinetic model 
incorporating the basic ingredients of non-ideal fluids with 
competing attractive/repulsive interactions. Although a 
one-to-one mapping with a corresponding physical sys- 
tem remains to be developed, the present model ex- 
hibits a number of highly non-trivial features of soft-glassy 
behaviour, including the Herschel-Bulkley rheology dis- 
cussed in this Letter. Finally, in light of the results dis- 
cussed in this paper, one could raise the following ques- 
tions: how far are present materials/experiments from 
the scenario depicted in this Letter? Can new mate- 
rials/conditions be adapted/designed in such a way as 
to realize the scenario revealed/suggested by the simu- 
lations? Since the present mesoscopic model can access 
scales close to experimental ones, we hope that the present 
work can raise new stimulating challenges for joint numer- 
ical/experimental work. 

* * * 

Valuable discussions with H.C. Oettinger, H.J. Her- 
rmann and I.V. Karlin are kindly acknowledged. 
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